#!/usr/bin/perl -w
use strict;

use FindBin;
use lib ("$FindBin::Bin/..", "/net/cpp-group/Leo/bin");
use db_parameters;
use ortho_species;
use run_cmd;

use DBI;


print STDERR <<"HEADLINE";
888888888888888888888888888888888888888888888888888888888888888888888

    table_micro_synteny

	Description:
		Obtain ENSEMBL IDs from the database or from the sequence file
	Populates DB tables:
		ens_id except for the PI field

888888888888888888888888888888888888888888888888888888888888888888888

HEADLINE

my $SQL_DISPERSED_FILTER = "reason_id & 4 = 0 AND ";
#my $SQL_DS_FILTER = "";
my $SQL_DS_FILTER = "o.dS < 1.09 AND ";


sub min
{
	my @t = sort {$a <=> $b} @_;
	return $t[0];
}

#-----------------------------------------------------------------------------------------
#
#	Read in $species micro synteny
#
#-----------------------------------------------------------------------------------------
sub load_micro_synteny($$$)
{
	my ($dbh, $protocol_id, $analysis_set) = @_;
	print STDERR "\tRead in micro synteny\n";
	
	$dbh->do("COPY orthologs.micro_syntenic_blocks(species, protocol_id, analysis_set, ".
			"micro_syntenic_block_id, disentangled_gene_id) FROM STDIN") or die DBI::errstr;
	for my $species($ortho_name1, $ortho_name2)
	{
		open (MICROS,  "$dir_orthologs_output/synteny/$analysis_set.micro.$species.list") or die;
		$dbh->pg_putline("$species\t$protocol_id\t$analysis_set\t$_") for (<MICROS>);
		open (ORPHANS,  "$dir_orthologs_output/synteny/$analysis_set.$species.orphans") or die;
		$dbh->pg_putline("$species\t$protocol_id\t$analysis_set\t$_") for (<ORPHANS>);
	}
	$dbh->pg_endcopy() or die DBI::errstr;
}

#-----------------------------------------------------------------------------------------
#
#	Read in $species macro synteny
#
#-----------------------------------------------------------------------------------------
sub load_macro_synteny($$$)
{
	my ($dbh, $protocol_id, $analysis_set) = @_;
	print STDERR "\tRead in macro synteny\n";
	
	$dbh->do("COPY orthologs.macro_syntenic_block_size_histograms(species, protocol_id, analysis_set, ".
				"block_size, cummulative_gene_fraction, count) FROM STDIN")  or die DBI::errstr;
	for my $species($ortho_name1, $ortho_name2)
	{
		open (MACROS,  "$dir_orthologs_output/synteny/$analysis_set.macro.$species.histogram") or die;
		$dbh->pg_putline("$species\t$protocol_id\t$analysis_set\t$_") for (<MACROS>);
	}
	$dbh->pg_endcopy()  or die DBI::errstr;
}






#-----------------------------------------------------------------------------------------
#
#	get_synteny_position
#
#-----------------------------------------------------------------------------------------
sub get_synteny_position($$$)
{
	my ($dbh, $protocol_id, $analysis_set) = @_;
	if ($analysis_set eq '1_to_1')
	{
		$analysis_set = '1 to 1'
	}
	else
	{
		$analysis_set = 'to'
	}
		
	my $sql_cmd =<<"PL/SQLCMD";
    SELECT 
        species, 
        chromosome, 
        ortholog_id, 
        prot_position, 
        CASE WHEN strand THEN '1' ELSE '0' END, 
        disentangled_gene_id 
    FROM 
        taxon.prot_positions NATURAL JOIN 
        orthologs.ortholog_sets o NATURAL JOIN 
        orthologs.gene_pseudo_evidence NATURAL JOIN 
        taxon.gene_loci
    WHERE 
        $SQL_DS_FILTER
        $SQL_DISPERSED_FILTER
        ortholog_type ~ '$analysis_set' AND 
        chromosome ~ '^[0-9XY]+\$' AND 
        protocol_id = $curr_protocol_id
    ORDER BY 
        prot_position
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd."\n";
	return $data
}



#-----------------------------------------------------------------------------------------
#
#	load_table_micro_syntenic_blocks
#
#		load the micro_syntenic_blocks table
#
#-----------------------------------------------------------------------------------------
sub load_syntenic_blocks_table($\%)
{
	my ($dbh, $protocol_id_to_name) = @_;
	
	$dbh->do("DELETE FROM orthologs.micro_syntenic_blocks WHERE protocol_id = $curr_protocol_id;");
	$dbh->do("VACUUM ANALYSE orthologs.micro_syntenic_blocks;");
	$dbh->do("DELETE FROM orthologs.macro_syntenic_block_size_histograms WHERE protocol_id = $curr_protocol_id;");
	$dbh->do("VACUUM ANALYSE orthologs.macro_syntenic_block_size_histograms;");
	my @protocols = get_orthologs_protocols();
	push(@protocols, ['1_to_1', 0]);
	for my $protocol (@protocols)
	{
		my ($protocol_name,  $protocol_id) = @$protocol;
		next unless $protocol_id == $curr_protocol_id;
		$protocol_id_to_name->{$protocol_id} = $protocol_name;
		for my $analysis_set('1_to_1', 'all')
		{
			print STDERR "\tOrthology protocol [$protocol_name]: $analysis_set...\n";
			{
				print STDERR "\t\tGet syntenic loci for [$protocol_name]\n";
				my $data = get_synteny_position($dbh, $protocol_id, $analysis_set);
				open (FILE, ">$dir_pipeline_temp/synteny.data") or die;
				for (@$data)
				{
					print FILE join ("\t", @$_), "\n";
				}
			}

			print STDERR "\tBuild synteny for [$protocol_name]: $analysis_set\n";
			my $cmd1 =<<"bourneshellCMD";
			$dir_lg_share_bin/build_syntenic_blocks
				--remove_singletons
				--orphans $dir_orthologs_output/synteny/$analysis_set.[SPECIES].orphans
				--macro_synteny_histogram $dir_orthologs_output/synteny/$analysis_set.macro.[SPECIES].histogram
				--micro_synteny_histogram $dir_orthologs_output/synteny/$analysis_set.micro.[SPECIES].histogram
				--list_micro_synteny $dir_orthologs_output/synteny/$analysis_set.micro.[SPECIES].list
				--chromosome_macro_histogram $dir_orthologs_output/synteny/$analysis_set.macro_per_chromosome.[SPECIES].histogram
				--full_macro_synteny $dir_orthologs_output/synteny/$analysis_set.full.[SPECIES].macro.synteny
				--full_micro_synteny $dir_orthologs_output/synteny/$analysis_set.full.[SPECIES].micro.synteny
				--verbose  < $dir_pipeline_temp/synteny.data
bourneshellCMD
			run_cmd($cmd1);

			load_micro_synteny($dbh, $protocol_id, $analysis_set);
			load_macro_synteny($dbh, $protocol_id, $analysis_set);
		}
	}
	$dbh->do("VACUUM FULL orthologs.micro_syntenic_blocks;");
	unlink ("$dir_pipeline_temp/synteny.data");
}




#-----------------------------------------------------------------------------------------
#
#	load_table_micro_syntenic_blocks
#
#		load the micro_syntenic_blocks table
#
#-----------------------------------------------------------------------------------------
sub load_table_micro_syntenic_block_size_histograms($)
{
	my ($dbh) = @_;
	#
	#	load block sizes
	#
	$dbh->do("DELETE FROM orthologs.micro_syntenic_block_sizes WHERE protocol_id = $curr_protocol_id;");
	$dbh->do("VACUUM ANALYSE orthologs.micro_syntenic_block_sizes;");
	$dbh->do("DELETE FROM orthologs.micro_syntenic_block_size_histograms WHERE protocol_id = $curr_protocol_id;");
	$dbh->do("VACUUM ANALYSE orthologs.micro_syntenic_block_size_histograms;");
	$dbh->do("CREATE TEMP TABLE t1 ( block_size INT4, count INT4, species TEXT, protocol_id OID, analysis_set TEXT)");
		my $sql_cmd =<<"PL/SQLCMD";
        INSERT INTO orthologs.micro_syntenic_block_sizes 
            (block_size, species, protocol_id, micro_syntenic_block_id, analysis_set)
            SELECT 
                    COUNT(*), species, protocol_id, micro_syntenic_block_id, analysis_set
                FROM 
                    orthologs.micro_syntenic_blocks 
                WHERE
                    protocol_id = $curr_protocol_id
                GROUP BY 
                    protocol_id, species, micro_syntenic_block_id, analysis_set;
        INSERT INTO t1
            ( block_size, count, species, analysis_set )
            SELECT 
                    block_size, count(*), species, analysis_set
                FROM
                    orthologs.micro_syntenic_block_sizes
                WHERE
                     protocol_id = $curr_protocol_id
                GROUP BY 
                    species, block_size, analysis_set;
PL/SQLCMD
	$dbh->do($sql_cmd) or die DBI::errstr.$sql_cmd."\n";


	#
	#	Calculate and load cummulative block sizes
	#
	$sql_cmd =<<"PL/SQLCMD";
      INSERT INTO
            orthologs.micro_syntenic_block_size_histograms
            (
            block_size, 
            count,
            species, 
            protocol_id,
            cummulative_count, 
            analysis_set
            )
        SELECT 
                b.block_size, 
                b.count, 
                b.species, 
                $curr_protocol_id, 
                SUM(a.block_size * a.count) AS cummulative_count,
                b.analysis_set
            FROM  
                t1 a, 
                t1 b 
            WHERE 
                a.block_size <= b.block_size AND 
                a.species = b.species AND
                a.analysis_set = b.analysis_set
            GROUP BY 
                b.species, 
                b.block_size, 
                b.count,
                b.analysis_set
            ORDER BY 
                b.analysis_set, 
                b.species;
PL/SQLCMD
	$dbh->do($sql_cmd) or die DBI::errstr.$sql_cmd."\n";
	$dbh->do("VACUUM FULL orthologs.micro_syntenic_block_size_histograms;");
}


#-----------------------------------------------------------------------------------------
#
#	load_table_micro_syntenic_block_size_summary
#
#		load the micro_syntenic_block_size_summary table containing N50 sizes for the different 
#		protocols
#
#-----------------------------------------------------------------------------------------
sub load_table_micro_syntenic_block_size_summary($\%)
{
	my ($dbh, $protocol_id_to_name) = @_;
	#
	#	Get medians
	#
	my $sql_cmd =<<"PL/SQLCMD";
	SELECT 
	    MIN (block_size), 
        a.species,
        a.analysis_set
    FROM 
        orthologs.micro_syntenic_block_size_histograms a, 
        (
            SELECT 
                    species, 
                    MAX (cummulative_count / 2) AS median,
                    analysis_set
                FROM 
                    orthologs.micro_syntenic_block_size_histograms 
                WHERE
                    protocol_id = $curr_protocol_id
                GROUP BY 
                    species,
                    analysis_set
                ORDER BY 
                    species,
                    analysis_set
        ) AS b
        WHERE 
            a.cummulative_count >= b.median AND 
            a.protocol_id = $curr_protocol_id AND 
            a.species = b.species AND
            a.analysis_set = b.analysis_set
        GROUP BY 
            a.species,
            a.analysis_set ;
PL/SQLCMD
	my $data = $dbh->selectall_arrayref($sql_cmd) or die DBI::errstr.$sql_cmd."\n";

	# save to panda
	$dbh->do("DELETE FROM orthologs.micro_syntenic_block_size_summary WHERE protocol_id = $curr_protocol_id;");
	$dbh->do("VACUUM ANALYSE orthologs.micro_syntenic_block_size_summary;");
	$dbh->do("COPY orthologs.micro_syntenic_block_size_summary(protocol_id, protocol_name, median_block_size, species, analysis_set) FROM STDIN");
	my $protocol_name = $protocol_id_to_name->{$curr_protocol_id};
	
	$dbh->pg_putline(join ("\t", $curr_protocol_id, $protocol_name, @$_)."\n") for (@$data);
	$dbh->pg_endcopy() or die DBI::errstr().".";

	# print as well
	print "\nN50 block sizes:\n";
	for (@$data)
	{
		print ("\t". join ("\t", @$_). "\n");
	}
	print "\n\n";
};





#
# connect to genome database
#
my $dbh = connect_to_panda();
use get_psql_time;
my $start_time = get_psql_time($dbh);
$dbh->{RaiseError} = 0;

system ("mkdir -p $dir_orthologs_output/synteny");

#
# load tables
#
my %protocol_id_to_name;
load_syntenic_blocks_table($dbh, %protocol_id_to_name);
load_table_micro_syntenic_block_size_histograms($dbh);
load_table_micro_syntenic_block_size_summary($dbh, %protocol_id_to_name);

print STDERR "\tCompleted\n";



log_pipeline_stage($dbh, $curr_protocol_id, 42, 'table micro/macro synteny_blocks', $start_time);
$dbh->disconnect();





